//
// Created by sunmy on 2020/10/16.
//

#ifndef HW1_GAUSS_ELIMINATION_H
#define HW1_GAUSS_ELIMINATION_H

#include "Eigen/Core"
#include <iostream>

template <int n = 10, typename T = long double>
void gauss_elimination_impl()
{
    using namespace Eigen;
    using namespace std;
    typedef Matrix<T, n, n> MatrixND;
    typedef Matrix<T, n, 1> VectorND;

    MatrixND H = MatrixND();

    for (auto i = 0; i < n; ++i)
    {
        for (auto j = 0; j < n; ++j)
        {
            H(i, j) = (T)1.0 / (i + j + 1);
        }
    }

    VectorND b = H * VectorND::Ones(n, 1);

    //    cout << "Now the matrix H is:" << endl << H << endl;
    //    cout << "Now the vector b is:" << endl << b << endl;

    for (auto i = 0; i < n; ++i)
    {
        auto maxarg = i;
        auto maxvalue = abs(H(i, i));
        for (auto j = i; j < n; ++j)
        {
            if (abs(H(j, i)) > maxvalue)
            {
                maxvalue = abs(H(j, i));
                maxarg = j;
            }
        }
        if (maxarg != i)
        {
            // swap two lines
            H.row(i).swap(H.row(maxarg));
            T t = b(i);
            b(i) = b(maxarg);
            b(maxarg) = t;
        }

        // regularize line i
        for (auto k = i + 1; k < n; ++k)
        {
            H(i, k) = H(i, k) / H(i, i);
        }
        b(i) /= H(i, i);
        H(i, i) = 1;

        for (auto j = i + 1; j < n; ++j)
        {
            for (auto k = i + 1; k < n; ++k)
            {
                H(j, k) = H(j, k) - H(i, k) * H(j, i);
            }
            b(j) -= b(i) * H(j, i);
            H(j, i) = 0;
        }
    }

    VectorND y = VectorND();
    y(n - 1) = b(n - 1);
    for (auto i = n - 2; i >= 0; --i)
    {
        y(i) = b(i);
        for (auto j = i + 1; j < n; ++j)
            y(i) -= y(j) * H(i, j);
    }

    cout.precision(12);
    cout << "The solution for n = " << n << " is:" << endl
         << y << endl;
}

void gauss_elimination()
{
    gauss_elimination_impl<10>();
    gauss_elimination_impl<11>();
    gauss_elimination_impl<12>();
    gauss_elimination_impl<13>();
    gauss_elimination_impl<14>();
    gauss_elimination_impl<15>();
}

#endif //HW1_GAUSS_ELIMINATION_H
